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ABSTRACT 

This work deals with the dynamical stability of flexible 
cylindrical rods subjected to external parallel, uniform and steady 
flow* A general finite element model is presented for isolated rods 
considering old and new formulations. System matrices are obtained 
using finite element based Galerkin ' s method and natural boundary 
conditions have been considered explicitly in these matrices along 
with the geometric boundary conditions* The results obtained in the 
form of oomplex eigenvalues for old and new formulations have been 
compared, wherever possible, with the available results through 
Argand diagrams, and discussed* Good agreement is found for pinned- 
pinned and cantilever (with tapered free end) rods except in the ill- 
conditioned range, where even available results can not be relied 
upon* Results for clamped-c lamped rods for old f ozmulati cn are 
corrected* New results are shown for clamped-c lamped rods for new 
formulation* 

A generalized computer programme is developed to form the 
dynamical matrix, vfoich gives complex eigenvalues by calling a library 
subroutine 2 INDIA. An important feature of this programme is its 
flexibility. It can accommodate all the boundary conditions by 
making few changes in the input cards. 



CHAPTER I 


INTRODUCTION 

1.1 INTRODUCTION TO RIOT-INDUCED VIBRATION PRQBIBM 

The response of a structure, when subjected to an external 
flow and/or containing a flowing fluid, is of considerable practical 
importance and is encountered in various engineering system components, 
such as nuclear fuel bundles (rods) and heat exchanger tubes etc. 

The dynamic behavior of such structures is quite complex and it has 
been a matter of interest for many dynamicists for last two decades. 
This type of problem, where the flow of the fluid, external and/or 
internal, is responsible for the vibration of the structures, is 
termed as "Plow-Induced Vibration Problem". In reference to the 
orientation of the fluid flow with respect to the axis of the struc- 
ture , the flow induced vibration problem may be classified as follows: 

1. Parallel flow, i.e., fluid is flowing parallel to the axis 
of the rod. It is further subdivided as: 

(a) Internal flow, 

(b) External flow, and 

(c) Internal and External flow. 

2. Cross flow, i.e., fluid is flowing perpendicular to the 


axis of the rod 
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Present work involves the study of dynamic stability of a 
single rod subjected to parallel external flow by the finite element 
method. 

Fluid flowing past a system is a source of energy, which is 
imparted to the system components. This energy is responsible for the 
vibration of the system components which may cause instability in the 
system. This instability, ultimately, may cause the system damage# 

At small flow velocities, instability usually does not occur* 
As the flow velocity increases, the structure becomes unstable at a 
oertain velocity# This velocity is called critical velocity* This 
instability may be of buckling type and/or oscillatory type. These 
instabilities are pronounced enough to cause failure of the system 
components, and hence need to be investigated so that the system can 
be designed to withstand the flow-induced vibrations. 

Analysis of the flow-induced vibration problems need the 
response, i#e., the actual motion of the structure, and the stability 
of the structure, ±#e*, to know, how sensitive the structure is for 
the changes made in the system parameters involved therein# The sta- 
bility analysis often supplies sufficient information to consider the 
vibration problem as solved. Moreover this analysis is comparatively 
simpler* The response analysis requires the complete knowledge of 
the fluid dynamics, i.ek., damping mechanism and flow originated 
excitation which may be due to unsteady pressure field or pulsations 


in the flow rate* 



3 


1.2 PREVIOUS W ORK ON THE DYNAMICS OE RODS SUBJECTED 
TO PAEAHEL EXTERNAL EIOW 


It was Burgreen et. al j_2 ] who first demonstrated the pheno- 
menon of rod vibration caused by the parallel flow of a fluid past a 
cylindrical rod in the year 1958, while studying the heat transfer of 
water flowing across a bundle of rods in a nuclear reactor core. He 
observed the oscillations to occur at all flow velocities. Also the 
oscillations were found to be relatively constant throughout the whole 
range of velocities which led him to believe that the vibration was 
of self -excited type. He developed an empirical expression for maxi- 
mum amplitude of vibration by dimensional analysis of differential 
equation of motion assuming excitation force to be proportional to the 
dynamic pressure p V ) of the fluid. This was followed by several 
experimental studies | 12, 15, 21, 22^J ? with different excitation 


forces, wherein the main aim was to study the character of the vibra- 
tions and obtain an empirical expression for maximum amplitude of vibra- 
tion which can be used to solve mechanical design problems in reactors 
and heat exchangers. 


In these experimental studies, various researchers assumed 
different forcing functions and obtained empirical expressions for 
mascimum amplitude of vibration. Most of them j~3 J selected turbulent- 
boundary layer pressure fluctuation as the forcing function. 

The first analytical and experimental treatment of the prob- 
lem was presented by Paidoussis fl4 1 in 1966. His mathematical 
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model was based on Bernoulli-Euler beam theory to describe the flexu- 
ral motions of the cylindrical body, slender body theory for the 
coupled inviscid hydrodynamic forces, and fairly reasonable relation- 
ships for the corresponding viscous forces. He considered three 
boundary conditions 5 namely pinned-pinned ends, clamped-c lamped ends 
and clamped-free ends -with tapered, streamlined downstream end and 
showed that small flow velocities damp the free motion of the system 
but high flow velocities induce buckling and oscillatory hydroelastic 
instabilities. He showed that flutter instabilities are caused by 
lateral frictional forces, and in the absence of hydrodynamic -drag 
effects only buckling instability is possible. Externally applied 
tension was observed to stabilize the system while compression des- 
tabilized it* 

Paid ous sis extended his previous study of the dynamics 

of flexible cylinders in axial flow and presented a general theory to 
account far the small, free, lateral motions of a flexible, slender, 
cylinderical body with tapered ends, totally submerged in liquid and 
towed at steady speed V# for particular shapes of the ends and length 
of ttw-rope , he showed that the body may be subjected to oscillatory 
and non-oscillatory instabilities for V > 0 5 at small V , these 
instabilities correspond to those of a rigid body. He also showed 
that at higher V, the system generally regains stability in the 
above modes, but may be subjected to higher— mode, flexural oscillatory 
instabilities* The effect of a number of dimensionless parameter an 
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stability were calculated to obtain the critical conditions of sta- 
bility. It was shown that optimum stability is achieved with a 
streamlined nose s a blunt tail and a short tow-rope. 

In 1973 , Paidoussis p7^j extended the theory by including 
material damping ? gravity, pressure forces and fluid damping force, and 
also applied the theory to cylinders in confined flows. In 'this work 
frictional forces were taken into account correctly. He showed that 
small flow velocities damp free motions , but at sufficiently high 
velocities hydroelastic instabilities occur. Pinned-pinned cylinders 
are subjected to buckling in their first mode , and then at higher 
flow velocities in their second mode? at slightly higher flow velo- 
cities coupled-mode flutter was observed, Critical velocities for 
the cylinder in a closely spaced cluster were much lower than those 
for an isolated one. Cantilever cylinders, on the other hand, were 
subjected to buckling in their first mode and to flutter in their 
higher modes. He carried out the calculations for cantilevers for 
high velocities, These enabled him to observe that the loss of sta- 
bility in the second mode is preceded by regaining of stability in 
the first mode? then in the similar manner second mode regains sta- 
bility at a slightly higher velocity than is necessary for the onset 
of third mode flutter. He also discussed the case of cylinders in 
axial flow subjected to an arbitrary force field. 

Rec sntly in i/75, Paidoussis 20 analysed the stability 


S' 

I 


of s Lind 
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Recently in 1975? Paidoussis J 20 j analysed the stability 
of slender cylinders with pinned-pinned ends and clamped-free ends, 
subjected to axial flow with harmonically perturbed velocity. He 
shewed the possibility of instabilities for certain ranges of frequen- 
cies and amplitudes of the perturbations. The parameters which des- 
tabilize the system in steady flow, make the parametric instability 
regions wider. In case of simply-supported cylinders parametric ins- 
tabilities occur over specific ranges of flow velocities. In the 
case of cantilever cylinders the instabilities were found to occur 
in some of the modes of the system and the regions of these instabi- 
lities in strutt diagrams were quite complex, Morewer, the regions 
of instability in strutt diagrams became wider as the system approached 
the point of instability in steady flow. 

1.3 PRESENT WORK 

In the present work the flow-induced vibration problem, in 
whioh vibration of the rods is caused by external parallel flow, has 
been solved by the finite element method and the dynamical behavior 
of the rods with different boundary conditions has been studied. 

Old f ormulati cn , Paidoussis Jj4j| , and, the new fozmulatiai, Pai— 
doussis Pi 7 ”| , both, have been considered for studying the dynamics 
of rods with pinned-pinned ends , clamped— clamped ends and clamped— 
tapered free ends. The results obtained in the present work have 
been compared, wherever possible, with the available results for 
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old and new formulations. Also, the effect of various non-dimensional 
parameters on the stability of tho erods has been studied and compared 
with the available results. 

In the second chapter the equations of motion are obtained 
in the matrix fora by the finite element analysis using Galerkin’s 
method. The stiffness, damping and mass matrices for a typical finite 
element are obtained. After assembling these matrices for the whole 
domain, natural boundary conditions (as a result of finite element 
based Galerkin's method) and geometric boundary conditions are intro- 
duced. Finally, the system dynamical matrix is obtained by arranging ^ 
these matrices. Method of solution is also discussed in this chapter. 

Complex eigenvalues of the dynamical matrix are computed 
using a library subroutine for various boundary conditions. The 
various dimensionless parameters, whose effects studied on the sta- 
bility of pinned-pinned rods are s 

(a) effect of dimensionless mass ( g ) 

(b) effect of external tension (F) 

(c) effect of internal damping (a) 

(d) effect of gravity ( y) 

All the results obtained arc illustrated and discussed in 
the third chapter. Furthermore, mechanism of instability is also 
discussed briefly in this chapter. 
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Chapter four contains the conclusions on the basis of the 
results obtained in the present v/ark. 

Further, in Appendix A, various fluid forces resulting from 
the fluid flow have been discussed in detail. Also governing equa- 
tion for small lateral motions for an isolated rod has been derived. 

Appendix B deals with the old formulation of the problem 
and its finite element analysis. The difference between old and new 
formulations has been discussed in brief. Appendix C gives the boun- 
dary conditions for a cantilever rod with tapered free end. 

lastly, the computer programme has been introduced in 
Appendix 1 whose detailed listing is given in Appendix E» 



CHAPTER II 


R0RM7IATI0H OR PR OB IBM AND COMPUTATIONAL APPROACH 


In the present chapter the fifth order differential equa- 
tion governing the transverse motion of a rod in an infinite fluid 
flowing at a velocity V is converted into simultaneous differential 
equations by finite element method using G-alerkin's method, further, 
these equations are rearranged to form standard eigenvalue problem 
to obtain complex eigenvalues. 


2.1 FINITE ELEMENT ANALYSIS 5 

The governing equation of motion for the small transverse 
vibration of a flexible slender rod of circular cross-section immersed 
in an infinite fluid flowing with a uniform velocity V can be 
written as, Paidoussis | 17 J , 
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with the boundary conditions as 


El *Mr + ul *Mr ( I + (l-2v)(pA.))|^ 


+ pi — 2 

ax J ax ' 1 at 


+ K y = 0 or y = 0 at x = 0, 
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( 2 . 2 ) 


where Z and W are given by Bqns. (2.3) and (2.4), respectively^ 
and 


y s internal damping coefficient representing a Kelvin-Voigt 

type of damping in the material of the rod. 

El s flexural rigidity of the rod. 

C ss coefficient of virtual mass, 
m 

m^ s= virtual mass of rod per unit length. 

m = mass of rod per unit length, 

r 

V ss mean flow velocity of the fluid. 

6 ss. equal to 0 or 1 accordingly as the rod is free or net 
to slide axially at x = £• 
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dL 


to, 

tli 


externally applied, uniform tension. 

Poisson’s ratio. 

steady-state pressure of the fluid. 

cross-sectional area of rod* 

pA at x = * 

density of the fluid. 

diameter of rod. 

skin friction drag coefficient. 

coefficient of viscous damping. 

hase drag coefficient, 

length of rod. 

axial co-ordinate. 

time , 

lateral displacement of rod. 

displacement spring constants at x = 0 and x = L, respectively, 
rotational spring constants at x = 0 and x = L, respectively. 


Por the sake of simplicity, various terms in Eqn. (2.1 ) 

2 

associated with ■'Mj- , and ^ are put in the following manner : 

9 X 3 u 


3x 


Z o (C m m f V 2 -{ <S [ T + (1 - 2 v) ( pA ) ]■+[•“ p f P V 2 0 f + Gya^g] 
1(1 X*] + ( \ P f D 2 V 2 C b (1 - 5)]}) ? (2,3) 
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N. 

i 


*2 (e) .5 (e) .2 (e) 

(m + C m ) ^ ~ + y I ~ + 2 C m V 

r m f 3t 2 3x 4 3t " * 


"m f 3 x 3t 


M 


*4 (e) 


s? ■ ( e ) 


+ X 


3 t 


+ El 


+ Z 


3 x 


+ W ((x +x) 


3 2 (e) a (e)~ 

£_x_ . ix- ) 

2 + 3 x } 


3 x“ p 3x‘ 

i = 1 ,2, ,f. (2.8) 

Now, the above equation is integrated by parts to introduce 

the natural boundary conditions and reduce the order of derivatives 

s imultane ous ly . 


dx=0, 


Second term is integrated twice and it results, 


»5 (e) 


BT. V I - 

0 1 3x 4 3t 


dx = 


N. VI _ 
x _3 


4 (e) 

tlLL 


ax at 


33 (e) .h 

ax V1 rk~ 


d®. 

1 


r h !l H i 

+ J 

0 dx 


VI 


s 3 (e) 

ix—ax . 

3/3t 


3x 3 t I 0 
(2.9) 


Third term, after integrated once, yields, 


h g2 (e) 

f ®. 2 C m_ V dx = 

J Q l mf 3x 3 1 

h d®. 


3 >> 

If. 2 C m _ V m 
x m f 3t 


0 


dx 


2 C nn V 


m f 3 1 


9v ( e ) 

3-Z dx . 


(2.10) 


Fifth term is also integrated twice and it gives, 

a 3„(e) , h i dSL 


h |\4 (e) 

/ ®. El — ^~r — dx = 

0 1 3x 4 


®. El 5- 

i ^ 3 


3x' 


0 


3 2(e),h 

— - El 
dx - 2 


3 x 


0 


2 

h d ®. 


+ J 2 
0 dx 


i -ex <Lz 


2 (e) 


El 


3x 


dx . 


( 2 . 11 ) 
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Sixth term reduces to, 


h 3 2 v ( e ) 

f N. Z 

] 0 1 9 x 2 


dx 


IZ 


(e) 


I. Z 
1 9x 


-/ 


h an 


i Z 3_X 

dx 3 x 


(e) 


dx 


0 0 

( 2 . 12 ) 

(ttoo ieX.>».j osus. 7 Gc(_J 


Seventh term turns out to be very interesting^ and its 


integration by parts gives, 

f h , (e) (e) 

; 0 1, i w((x p + x) 3 -^~ + 3 i-" 

h an 


dx 


JT. W (x + x) 
x p ' 9x 


(e) 


h 


0 


/ ^ » (x. + x) 


§x 


(e) 


3x 


dx 


(2.13) 


Q v ~*‘ P 

Substituting Eqns. (2.9), (2.10), (2.11), (2.12) and (2.13) in Eqn. 
(2.8), it takes the following form, 


h 

J q {[ U r + C m m f ) n ± ' ] + [y I ~~~2 ^ 2 — - w m “ f 'cbc 


^ (e) 


d U. a 3 (e) 


an. 


rj v v. ,. hlT n _ 2 c mJ z±4z i 


(e) 


dx 9x 9t 


r (e) 


d ll. ,2 ( e ) 


an. 


+ X *i 9t 


B^z. ] + [EI i^i^. z ri ax 


(e) 


a 2 . 2 
dx ox 


dx 9x 


an. „ (e) 
r _xaxU. 

p dx & 


an. ^ (e) 

,} te + 


N. 2 C m V 
x m ± at 


a 4 (e) h 
S. p I S -X— ' 

1 3x 3 3t 

i£± 


an. 3 (e) ,h 
1 y 


dx 


9x 2 9 t 


0 


0 


n. ei 

1 3x 5 


an. a 2 (e) 

1 O V 


dx 


EI 


9 x 


-b 


n z 

i 9X 


(e) 


h 


0 


n ± w (x p + x)!j-X 


(e) 


h 


= 0 


X~ £ * 


(2.14) 



15 


Again, to have the system matrices, Eqn. (2.7) is substituted 
into the above E q n. (2.14) and rearranging, one gets for the element, 


fm] (e) {y} ne +[ c] (e) {£} Ue +[k] ( e ) {y} + 


1L 1114 
1 9- 3 


4 (e) 




a 4je) t ,h 
v y I 2- X- 


3x' 


0 


w i 


izfl 


3 x 9t 


h 


f a 3 T (e) ,h 

N 2 yl 

ax 3 1 


r > 


1L 2 0 m V 
1 m f 3 t 


(e) , If'- 


a 4 y ( e ) 

ox t 


h 

0 




,, 2 2 0 m ” f 7 li 


(e) , h 


0 


' 3 3 v ( e ) h 

1-T ux 

f M x „ 2 „ 

a x a to 


i 


2 C m V 
f m f a t 


(e) 


h 


“V T (e), hi 

I 1 El — X- - 




ax' 


0 




<\3 (e) I h 

K BI 

2 fe 3 1 0 


*>» 


X 


> a 2 J e ) h 

N 1 El 

ax 0 

2 , r (e) h. 

0 


n' ei -f 

3x 


a 3 (e)i h 
BI ~-Z — - 


ax 






r a 2 (e) j h 

XT EI 

8x 1 0 
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IL Z 


Z 


„ (e) 

axlx 

3x 

h ; 

i 

0 


1 

t | 
1 

"1 ® (x p + 

3yi e) 

h 

| 


S 2 W (x p + x) || 

3x 

0 




1 


s 


.(e) 


.(e) 




IX 


W f 2 gx 


(e)i 


h 

0 


1 

where 

Em] 


N f ¥ (x p + x) S 


(e) 




(e) 


Ec] 


(e) 


= mass matrix for the element 
h 

= / (m +G mj {1} I'S'I & 

0 r m f 

= total damping matrix for the element 

= Ecl] (e) +Ec 23 (e ) + Ec 3 l (e) 


M 


[cl]' 7 = internal damping matrix for the element 

h 


= / )i i {n } L N J ax 

0 


h 

0 

h 

0 

h 

0 


>=• (0} 


i(e) 


t c2 ] v ' s coriolic force dampirg matrix for the element 


la 

= ~ / 2 0 m V { N } { F J dx 

4 q m x 


i(e) 


[03]'^ = viscous drag damping matrix for the element 

h 

= / X { I } |_N J dx 

0 


(2.15) 


( 2 . 16 ) 


(2.17) 


(2. 18) 


( 2 . 19 ) 


( 2 . 20 ) 
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[k] 


(e) 


s= total stiffness matrix: for the element 


[kl] 


(e) 


:[ k1] (e) - [k2]^ - [k3] (e) 


elastic stiffness matrix for the element 


h 


( 2 . 21 ) 


=/ El{ I } (_N J dx 
0 


(2.22) 


( e \ 

[ k2] v 1 — stiffness matrix for the element due to curvature term 


a 2 

O y 

'■ mmi' jm ■ m 

8x 2 

h f , 

/ (Z + Wx){N}|jrj dx 
0 P 


(2.23) 


and 


[k3l^ = stiffness matrix for the element due to variable 

A 

curvature term x —fj 

3x 

h it 

= / fit I ) U J dx 
0 

( e ) 

To calculate these matrices one has to choose y . The 

(e) 

following conditions must be fulfilled while selecting this y 
to ensure monotonic convergence , Huebner f 9 ] : 

(i) Completeness , i.e., all uniform states of y and its partial 

derivatives upto the highest order appearing in Eqn. (2.14) 

(e) 

should have a representation in y J when in the limit, the 
siae of the element shrinks to aero. 


(2.24) 
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(ii) Compatibility, i.e 0 , at the interfaces of the elements y 


and any of its partial derivatives one order less than the 
highest order appearing in Kqn* (2.14) must be continuous. 

The following shape function, usually used for beam problems, 
in terms of generalized co-ordinates has been chosen for the present 
problem, since it satisfies the above requirements of completeness 
and compatibility, respectively, 

LhJ = L5 1 2 (3-2S 1 ) E, 2 s 2 h K 2 (3 - 2 5 2 ) s/h], 


.... (2.25) 


where = 1 


’2 “ h * 


Introducing the shape function into the Eqns . (2.16), (2*18), 
(2.19), (2.2o) , (2.22), (2.23) and (2.24) and integrating, one obtains 


[ cl] 


(m + C m )h 
r m r_ 

420 


(e) U I 


156 22h 54 -13h 

2 2 
4h I3h -3h 


Symmetric 


12 6h 


4h -6h 


Symmetric 


156 -22h 


(2.26) 


(2.27) 
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[c2] 


(e) 


2 0 m_ V 
xn f 

60 


[c3] 


(e) 


3k. 

420' 


t«l (e) = 3 

h 


(e) 


30 

6h 

30 

-6h 



-6h 

0 

6h 

-h 2 



-30 

~6h 

-30 

6h 


(2.28) 

6h 

h 2 

-6h 

0 





156 

22h 

54 

-13h 



2 

4h 

13h 

-3h 2 

(2.29) 


Symme trie 


t k2 = - 30h 


and. 


[k3l^ W_ 

= 30 

1 

Symmetric 

Therefore Eqn. (2 . 1 5 ) becomes . 


156 


-22h 

2 

4h 


3h 


12 6h - 

-12 

6h 

2 


2 

4h 

-6h 

2h 


12 

-6h 

Symmetric 


2 

4h 

36 3h 

-36 

3h 

2 

4h 

-3h 

-h 2 


36 

-3h 

2 

Symmetric 


4h 

18 3h 

-18 

0 

2 

h 

-3h 

2, 

-h/ 2 


(2.30) 


(2.31) 


(2.32) 
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- U I 


< 


[a] (e) { y>“ + [ 0 ] (e) <y}“ + [k ,(O {y}I >e + 


- y i 


0 

3 V e ) 
3x 2 




+ 

> < 


“ 2 G m V 3 *X 

El I o X 


+ 


+ y I 


~ EI 


4 (e) 
-j . L 

3fc 3 3 t 
0 

a*v (e > 


iP 




3 

3x p 

3 t 

0 


<3 ( e ) 

r — 

Ph 


3 x' 


*3 (e) 

+ y I — 

. 3x 3 1 

X, 

r 


a 2 ( e ) 

SI i-£-- 


+ 2 0 m : v|Xl 

m f 3 t - “ 


r *< 


0 


0 


r3 (e)| 
+ EI -~£_ fa 
3x 2 


> 


V_ 






< 


+ EI 


3x 

0 

3 V e ) 
3x 2 


- Z 




(e) 


3 x 


~\ 


q. 




< 


w (x p+ x) ax_ 


+ z 


six 

9x 


(e) 


f 

< 


v 


( e ) 

+ W (x +x) 

p 3x 1 


(2.33) 


Eor the whole domain of the problem, Eqn. (2.33) becomes 
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where 


[ = assembled mass matrix for the whole domain before applying 

the boundary conditions 

[ C]^ = assembled damping matrix for the whole domain before apply- 


ing the boundary conditions 


and 


[ K]gg = assembled stiffness matrix for the whole domain before 
applying the boundary conditions «, 

For computational ease, Eqn. (2.34) can be rearranged as 


follows: 


r 


-2 C bi Yy 
m f J1 

0 
0 


[M] BB (y} n + [C ] BB ( y} n 


+ 




+2 C m Yy 
m X 

0 


m+1 


f 

S 




- z e 






1 


0 


0 



0 


0 



0 




+ 



• 

0 


\ 

: 

0 


+ Z 6 



+ ¥ m h 9 | 



m+1 



m+1 

0 



0 






J 


> 


23 


= - <v 


r 


t 


'N 

~3 (e) a 4 (e) 





- (BiU_ + ) 



0 


9 x 3 x t 

1 







32(e) a 3 v ( e ) 


0 



-(ei — - + y 1 - i" — ) ■ 





ax ax at 

1 

0 

J 

• 

\ * 


^ < 
hi 

0 


« 



0 


0 



• 

• 


~3 (e) J (e) 



* 


(El 9 | + Hi ") 



0 


ax ax a t 

m+1 







a 2 (e) a 3 (e) 


0 


, 

(ei-?— + : ) 


l 


ax ax t 

m+1 


} 

> 


(2.35) 


In Bqn. (2.35)> the values - 2 V and + 2 Y are 

added up with the elements 0^ ^ and ^gm+l 2m+1 '^ ie assem kled 
damping matrix , respectively. It should also be noted that the co- 
efficients of 0 ^ and 9 are summed-up with the elements ^ 

and K„ . _ .of the assembled stiffness matrix, respectively. 

2m+1 ,2m+2 

After substituting the boundary conditions, Eqn. (2.35) can 


be written as 

[M]{ y} n + [C]{y} n + [K] {y} n = < 0 > , (2.36) 

where [ M ] , E C ] and [K] are the assembled mass, damping and 
stiffness matrices, respectively, for the whole domain of the problem, 
after applying the boundary conditions* 


Equation of motion (2.1) derived by Paidoussis E 1 7 1 is much 
more general than that derived earlier, Paidoussis [14 ] , as it 
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takes into consideration the gravity and pressurization effects. 

In addition, it eliminated the error of non-inclusion of the fric- 
tional term F.^ ( 3y/sx) in the earlier work. Earlier formulation 
[ 14 ] is referred as "Old Formulation" and the later formulation 
t 17 1 used in this chapter as "New Formulation" in the present 
work and are given in Appendices A and B« 

The author found it fascinating to solve the old formulation 
problem also by the finite element method and is given in Appendix B, 

2.2 COMPUTATIONAL APPROACH i 

The solution of these homogeneous set of differential equa- 
tions is obtained as described in Meirovitch f 11 3 . As per this 

reference the inertia, stiffness and damping matrices should nece- 
ssarily be symmetric to arrive at the solution. But it is not essen- 
tial for the matrices [ M 3 , [ K ] and [ C ] to be symmetric , refer 

Frazer, Duncan and Collar [6 3 . 

Equation of motion of forced vibration of a damped system 
is 

[M] { ‘y (t) } + [ C ] { y (t) } + [ K ] { y(t) }= {Q(t)} , 

.... (2.37) 

Equations (2.37) represent a set of coupled ordinary linear 
differential equations which is difficult to decouple unless matrix 
[ C 3 is a linear combination of matrices [ K ] and [K] , which can 
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occur In particular cases only. For a general case, equations (2.37) 
can be reduced to a system of first order differential equations as 
follows , Meirovitch [11] , 

—1 

Premultiply ing Eqn. (2.37) by [ M ] 


Cy (t) > +[ M l” 1 [ G 3 { y(t) } + [M] - 1 [K]{ y(t) } = [ M 3 ^ {Q(t) > 


Pet { q(t) } 


1 [ 


and {l r (t) }=J 


f{°} 


( 2 . 38 ) 


[ 2n x 1 ] 


jf Q(t) *iE2nx 1 3 

(2.39) 


This leads to a set of 2n first order ordinary differential equa- 
tions 


E M 1 ' { q(t) } + E K 1 {q(t) } = {y(t) > , <2.40) 

where 



[ 0] 

[ M] 


-[ M] [0] 

[ M3 * 



and [ K 3 = 



E M3 

[ c 3 


[0] { K] 


are real matrices of order 2n. 

From Eqn. (2.40) we obtain the homogeneous set of differen- 
tial equations in matrix form as 

E M 3 { q(t)} + E£ 3 { q(t) } « { 0 } (2.42) 

The solution of Eqn. (2.42) is obtained by letting 


{ 9 (t)> = { 0 } , (s.«) 

where represents a vector consisting of 2n constant elements* 

Eqn. (2*45), when introduced in Eqn. (2.42), leads to the eigenvalue 


problem 


(2.45) 


which can be rewritten in the form 

1 


where 


[D]{^}= ~ { J# } , 


D = - [ X ] ~ 1 [ M ] 


[M3 


[0] 


[0] 


[0] 


[ K [ 


-1 


[0] [Ml 
[Ml [ c 3 


[ I] 


[K ]"" 1 [ m 3 


[k 3 —1 tel 


[2h x 2n] 


( 2 . 46 ) 


plays the role of dynamical matrix. The multiplication of parti- 
tioned matrices is handled as if the submatrices were ordinary ele- 
ments. The identity matrix [IJ in Egn. (2.46) is of order n. 


The eigenvalues of the above dynamical matrix (2.46) will 
in general be complex. Real parts determine the stability of the 
system and the imaginary parts give the frequency of vibration. If 
the real part is negative-, the system will be stable and vice-versa. 
The system will have oscillatory instability (flutter) if real part 
is positive and imaginary part is non-zero f whereas if the imaginary 
part is zero, the system will undergo buckling instability'. 

The formulation presented in this chapter is general and 
proves to be versatile because the effect of various boundary 
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conditions (geometrical and natural) can be studied easily* Also^ 
the effect of various system parameters such as velocity^ axial ten— 
sion^ fluid pressure ? gravity etc* on the dynamics of rods can be 
studied with special interest* 

A computer programme has been prepared to compute the eigen- 
values using built-in subroutine Z INDIA* Listing of computer programme 
is given in Appendix®* Computer programme has been discussed in 


Appendix TX 



CHAPTER III 


RESULTS AND RISCUSSIOH 


The results obtained in the form of dimensionless parameters 
are included and discussed in this chapter. Various boundary condi- 
tions of the rod and the effect of various system parameters on the 
stability of the rod are illustrated in detail. Complex frequencies 
of the lowest four modes have been obtained as a function of the 
dimensionless flow velocity and Argand diagrams have been plotted. 
All the results are given in terms of the dimensionless parameters. 
Lastly, the mechanism of instability has been discussed in brief. 


3 . 1 DliMSIOUIESS PARAMETERS 

following dimensionless parameters have been used in this 


work, 

3 = 

u = 

F = 

a = 


Y = 
n 


Dimensionless mass = m f /(m f + m^) 

1 /2 

Dimensionless flow velocity = (m^/El) VL 

2 

Dimensionless axial tension = T L /El 


Dimensionless internal damping coefficient 
= {1/ [ E (m + m ) ] } 1/2 y / L 2 


Dimensionless gravity parameter = (m^ - m ^) S 1 Z® 1 

. . 2 n 

Dimensionless pressure parameter = { P& ) ^ 
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Re(oj ) = Real part of the dimensi earless frequency ={ (m^ + m J/EX} 1 / 2 ^! 2 
Im( w) = Imaginary part of the dimensionless frequency 
= { (m f + m^/BI } 1//2 WjL 2 

a ^ = Dimensionless displacement spring factor = K I /El 

2 

a^. = Dimensionless rotational spring factor = K^L /El 

-J /P 

c^ = Dimensionless viscous drag coefficient = C^D(m^/Bl) ' 
c^ = Friction drag coefficient = (4/ fr ) C 
c^ = Base drag coefficient = ( 4 / ^ 

3.2 COWER&EHCE OF THE E IGENV AIDES ; 

Finite element analysis says that as the number of complete 
and compatible finite elements is increased, the results converge 
monotonically to the exact values. However, one has to choose finite 
number of finite elements. For the present work, first for every 
boundary condition results were obtained for a few non-dimensional 
velocities for various numbers of the finite elements. It was obser- 
ved that for five elements results were quite acceptable for all the 
cases. Sample results for pinned-pinned rod and cantilever rod are 
shown in Table 1 and 2,. for RED frcox .2 to- 7 and frem 2 to 5.* respec- 
tively , (pp .30 , 3.1 ) Jl t should be noted that the change in the frequencies 
after EEL = 4 is negligible. Therefore, it was decided to use five 
finite elements for determining the results. Moreover, this saved 


computational time also 














TABLE 
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3 . 3 OOMPIEX PRBQUEICIES OF R CTOS ; 

The dimensionless complex frequencies obtained for the pre- 
sent problem are plotted in Argand diagrams for different boundary 
conditions of the rod® The imaginary part of the dimensionless 
value Im( w ), which represents the frequency of vibration is plotted 
on the abscissa and the real part He (to ) 5 on the ordinate* The various 
parameters involved in the problem are ct, 3 , Y , F , , n , u, , a ^ , 

C £’ °b ^ V 

3.3.1 HmBD-mmSD ROB ; 

The boundary conditions for pinned-pinned rod are 

y = 0 at x = 0 and x = L, (3*1 ) 

r. 2 r, 3 

El 2-Z +y I — - o at x = 0, (3*2) 

3x 3x~ 3t 

2 3 

El + Ul = 0 at x = L. (3.3) 

3x 3 x 3 t 

Since pinned-pinned rods are of considerable practical inte- 
rest as compared to the cantilevered rods and clamped-clamped rods 5 
this boundary condition is studied in detail and the effect of various 
dimensionless parameters on the stability of the rod is discussed* 

3 c 3 • 1 * 1 O U) EOHHQMT I Ohi 

In this old formulation, transverse component of tangential 
frictional force is neglected, see Paidoussis [ H 3 and Appen- 

dix B. Pig. 3*1 shows the dimensionless complex frequencies of the 
lowest three modes of a system as a function of the dimensionless 
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flow velocity u , for- 6 = 0 = 1 , £ c^ = 1, T = 0, = 0^=1 5. These results 
were compared with the results of Paidoussis [14] and it was observed 
that results were matching with those of available results except for 
certain velocity ranges, i.e*, in the first mode for u > 3 and and 
in second mode for 5< u< 6« Careful study of the available results 
indicated that there was ill-conditioning for the above mentioned 
velocity ranges, i.e., for small changes in u, a large change in 
complex frequencies was observed. This led the author to believe 
that there is ill-conditioning in the system in these velocity ranges 
and the convergence is incomplete for five elements. Therefore, the 
results for these velocity ranges, where ill-conditioning was believed 
to be present, were studied by increasing the number of elements from 
2 to 8. These results are tabulated in Table 3o 

Still for 8 elements there is some difference in the present 
and Paidoussis 1 s [ 14 3 results in these velocity ranges. Moreover, 
the results of Paidoussis [14] may also not be reliable for the 
ill-conditioned ranges. 

3 . 3 . 1 . 2 M EW ECEMCJIi TI QCT : 

Case (A ) : 3 = 0.1 , e c f = 1 ,5 = = 1 , c d = a = Y = n = I = 0. 

Eig. 3.2 shows the dimensionless complex frequencies of the 
lowest three modes of the system as a function of the dimensionless 
flow velocity u. The results match very well with these obtained by 
Paidoussis [ 17 ] except for certain ranges of u, where ill-condition- 
ing was observed to be present, i.e., 6 < u< 6.5 in the second mode 
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and 7 < u <. 10 in the third mode. Studyirg Fig. 3,2 end Table 3 it 
was thought that not much will be gained by increasing the number of 
finite elements* 

It is observed that small flow velocities damp free motions 
of the rod* It is also noted that the behavior of the pinne d-pinned 
rod is changed a lot by the corrections made in the frictional force 
terms* In the older formulati on , Fig* 3*1 , buckling occurs in the 
first mode only and flutter occurs in second and third modes* Vhereas, 
in the new formulation, Fig* 3*2, buckling can occur in all the three 
modes and there is coupled mode flutter at u=; 6.48 (first mode for 
u > 6*48* It may be noted that there is no coupled-mode flutter for 
u > 6*48, as wrongly interpreted by Paidoussis [ 17 ] • 

Case (b) : $ = 0.48, e c f = 0.25 ? 6= 0^ = 1 , c^ = « = y= p=r = 0. 

T'ig* 3.3 shows the dimensionless complex frequencies of the 
lowest three modes for the system as a function of the dimensionless 
flow velocity u. Considerably good agreement of the results with 
those obtained by Paidoussis C 17 1 is observed. In the first mode 
ill-conditioning is present for 3< u<3.25 and 6.25< u< 6*4. Slight 
discrepency is observed in the results for 6.25<u <6.4. 

Results for $=0.1 and keeping the other parameters same 
are also shown in Fig. 3.3 for u upto 3.1. It is interesting to 
note that the increase in $ tends to increase the damping in the 
system, for small values of u. In fact, it can be seen from the 
coefficient of from Eqn. (2.6). 

d X 
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In this case second mode does not cross Im( (0 ) - axis. 

Thus no buckling occurs in the second mode. The first mode branch 
first crosses Im( to) - axis through the origin at u ^ 3 . 1 4 and again 
at u 2 6.2§. It leaves the real axis at us 6.295 and flutter occurs 
for u > 6,295. There is no coupled-mode flutter, as wrongly inter- 
preted by Paidoussis [ 17 ] « 

Case (0) ; EFFECT OF TENSION : g = 0.1 , e c = 0,25, r = 10, 

6 = C =1, C - a =Y = n=0. 
m 1 d 1 

Fig. 3.4 shows the dimensionless complex; frequencies of the 
lowest three modes for the system as a function of dimensionless flow 
velocity u. For this case, Paidoussis C 1 7 3 obtained the values of 
first buckling velocity, second buckling velocity and the critical 
velocity for coupled-mode flutter as u^ = 4.455, u^ = 7*04 and 
u q q = 7.11. The author obtained these values as 4.456, u^s 7*045 

and u 7.12, This slight discrepancy between the results is due to 
ill-conditioning in the system near these- velocities. It is observed 
that externally applied axial tension has stabilizing effect, i.e., it 
raises the critical flow velocities. 

Case (P) ; EFFECT OF INTERNAL DAMPING - ; g = 0.1, e c f = 0.25, 
a - 0.003 , ^ = = ^ ^ Y= *¥ -r= 0 . 

Fig<3.5 shows the dimensionless complex frequencies for the 
system as a function of the dimensionless flow velocity u. The 
dynamic behavior is very similar to the Case (A), see Fig. 3.2. 
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'The first and second buckling velocities, as obtained by Paidoussis 
[ 173 are = 3.143 and = 6.28 and the critical velocity 
for coupled-mode flutter is u = 6.52. In the present case, the 
author obtained these velocities as u^-j - 3.141, u^ 0 ~ 6.29 and 

u c qM 6.35. Slight discrepency in these results is due to ill-condi- 
tioning near these velocities. It is noted from Pig. 3.5 that coupled- 
mode flutter occurs at u 2 6.35 and flutter occurs in the first 
mode for u > 6.35. 

Case (E) : EEEBCT Cig GRAVITY : p = 0.1, sc = 0.25, y= 10, 

4 = C m = ’< °d =a=I l=r = 0* 

Pig. 3.6 shows the dimensionless corapler frequencies for the 
system as the function of the dimensionless flow velocity u« The 
first and second buckling velocities obtained by Paidoussis 1 17 3 
are = 3.107, = 6.30 and the critical velocity for coupled- 

mode flutter is u A = 6.44. The author obtained these values as 
cO 


u^i = 3.2, __ ~ c0 

is present for the velocity ranges 3*1 < u < 3.5 and 6 < u <6.5. 
It is observed that the inclusion of gravity parameter slightly 
reduces the values of Im(w ) while the effect on Re ( w ) is negli- 


— 6.3 and u _ a 6.45. Here also ill-conditioning 


gible 
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3.3.2 C lAIgSD— G LfllvIBSX) ROD : 

The boundary condition for elamped-elamped rod are 
y = 0 at x = 0 and x = L, (3.4) 

y’= 0 at x = 0 and x = L, (3.4) 

3.3. 2.1 O LD FORMimTIOlT : 

Fig. 3.7 shows the dimensionless c cmplex frequencies of the 
lowest three modes for the system as a function of the dimensionless 
flow Telocity u, for 8 = 0.1, e = 1 ) ^ = C^ = 1 , T - 0. Results 
obtained by Paidoussis [14] are also shown in the same figure. It 
is observed that the results do not match. In this case no ill- 
conditioning is observed in the present work or Paidoussis 's work 
[ 14 ] * Author’s results are correct as the deck used for pinned— 

pinned rod, modified for the clamped-clamped case was used. Correct 
results for u = 0 confirmed that the boundary conditions were 
applied correctly. 

It is observed that buckling will occur at u - 6*25 in "^e 
first mode which contradicts the earlier result, Paidoussis [ 14] • 

Buckling does not occur in second and third modes. Flutter Instabi- 
lity occurs in second mode at u - 7.5 end in third mode at u— 10.7* 

It is in terest ing to note that third mode regions stability at u* 11.6. 

3.3. 2. 2 HEW FORMJIATIOI ; 

Fig. 3.8 shows the dimensionless complex frequencies of the 
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lowest three modes for the system as a function, of the dimensionless 

flow velocity u, for g = 0.1,e cr, = 1 ,6 = C = 1, c, =a=Y = n = r 

f m d 

= 0. Comparing this with Pig. 3.7 it is seen that due to inclusion 
of transverse component of frictional force , P^ , the behavior 
changes considerably, flow instability occurs by buckling in its 
first mode at u ^6.25 and in the third mode at u « 12. And it is 
interesting to note that there is no flutter instability when new 
formulation is used. 

3.3.3 CANTIIBVBR ROD WITH f AiEBED PREE Eli t 

She boundary conditions for a cantilever rod with tapered 


free end 

are, see Appendix 0, 



y = y* = 0 at x = 0 

(3.6) 


3x 3 !i 5 !t f 3t ** 



3 2 

- (m r + f a f ) x e = 0 

( 3.7) 

soxd 

pi 2 i3 

BI + l* X if— - Q, an-I, 

(3.8) 


3x 3x 3t 


where 

1 

x% = *r J A(x) dx. 

1-1 


3. 3.3.1 

OLD PORMJ RATION: 



Pig. 3.9 shows the dimensionless complex frequencies of the 
lowest four modes for the system as a functions of dimensionless flow 



velocity u, for g= 0.5, £C f" 6 = 0, = 1, f =s 0.8, x e /L = 0.01, 


r= 0* Results match completely with those obtained by classical 
methods 5 Paidcussis [14 ] » Both the results are shown in Pig. 3.9 
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Here? it is important to note that the column matrices 


{ 

and{ 


H ± 2 m f T J— 
3 (e) 


S. ( EX hi 

1 sP 


} s { H Z 
4 (e) 




(e) 


3X 


U 


0 


H. 

1 


(x +x) _ 

p X 


Xe) 


h 


0 


+ PI 


I2L 


3 

8 * t 


} found in Bqn. (2.14) must 


x=L 


be taken into consideration in mass, damping and stiffness matrices, 
as shown by Deb l 26 3 . In the analysis these can not be neglected 
as suggested by Szabo and lee [ 24] . 


3 . 3.3.2 HEW PQRMJIATIOH : 

Pig. 3.10 shows the dimensionless complex frequencies of the 
lowest four modes for the system as a function of dimensionless flow_ 
velocity u, for 0 = 0 . 5 , e c f = 1 , 5 = 0 , C^ = 1 , f = 0 . 8 , x e /L = 0 . 01 , 
= a=yri= Y= c^ = 0. Complete agreement was observed with the results 
of Paidoussis [171 upto u = 8 . After u> 8 , ill— conditioning was 
observed and for this velocity range, the convergence was studied by 
increasing the number of finite elements from 2 to 9. Results are 
tabulated in Table 4. Also in Pig. 3.10, for u > 8 , results with 
HE I _ 5 , HE1 - 8 and the results obtained by Paidoussis [ 17] are 

shown. Satisfactory convergence is not observed for u> 8 because 
of ill-conditioning. 

Comparing the Pigs. 3.9 and 3.10, it is observed that there 
the behavior of the system by including the 


is no basic change in 
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transverse component of the tangential frictional force in the 

new formulation* Buckling occurs in first mode at u^ =j 2.04 and 
flutter occurs at u^ - 5»15 in the second mode and at « 8,25 
in the third mode. 

It should be noted that the results can not be relied upon 
for velocities u > 8 , as ill-conditioning occurs in this velocity- 
range. 

3,4 THE MECHANISM OP INSQABIIITY ; 

After studying the various results , the mechanism of various 
hydroelastic instabilities for rods subjected to external flows may 
be explained in brief as follows s 
3,4,1 BUCKLING INSTABILITIES : 

2 g 8 

It can be seen that the term m„ If g— in Z g ^ ie 

1 S x 3x 

gover ning equation of motion, Eqn* (2,6), behaves more or less as an 

equivalent compressive force which increases with the increase in 

flow velocity? thus reducing natural frequencies of the system. It 

is possible that for certain velocity, frequencies may became zero 

causing buckling instability.- 

Eor p inn ed— pinned and clamped-clamped rods one may regard 
2 

the tena 2 4 as a generalized centrifugal force; when this force 
overcomes the flexural restoring forces, the rod buckles. Eor canti 
levered rods with tapered free end buckling may occur when the li ft 
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force generated at the inclined tapered free .end overcomes the restor- 
ing flexural forces. Although these effects are modified by the pre- 
sence of viscous forces, they remain dominant so far buckling insta- 
bilities are concerned, Paidoussis [ 1? ] . 

3.4.2 PUTTIER IHSTABIhlTIES s 

Oscillatory instabilities can be explained by considering 
the energy balance. When the energy transfer from fluid to rod and 
vice-versa, in the course of one cycle of oscillation, exactly balance, 
the condition is said to have neutral instability, i.e., dynamic 
equilibrium. But when the energy transfer from fluid to the rod is 
more than the energy transfer from rod to the fluid, the amplitude 
increases without limit causing flutter. In the opposite case osci- 
llations are damped. 

Paidoussis and Issid 27 have shown that conservative 
systems can only be subjected to coupled - mode flutter type of ins- 
tability, i.e., frequencies of two modes are equal, in addition to 
buckling type of instability. But author finds that conservative sys- 
tem can also be subjected to uncoupled flutter, see Pigs. 3.2, 3.3, 
3.4, 3.5 and 3.6. It is left for future workers to explain this.' 

Por detailed study of mechanism of instability one may see 
Refs. I 14, 17, 2? 3 . 


CHAPTER IV 


GOEOLUSIOHS 

In this chapter 9 the summary of the analysis made and the 
conclusions arrived at are illustrated as follows: 

1 • The finite element models developed for the new and old 

formulations are correct as the results obtained for the above 
mentioned formulations match quite well with those obtained 
earlier by classical methods for the cases of pinne d -pinne d rods 
and cantilever rods with tapered free end. 

However, for certain ranges of the flow velocities , ill- 
conditicning is present in the system for the pinned-pinned rods 
and cantilever rods with tapered free end. Eor these velocity 
ranges, the results obtained by finite element method and cla- 
ssical method may not be reliable . 

By increasing the number of finite elements the results do 
not converge for ill-conditioned velocity ranges for the cases of 
pinned-pinned rods and cantilever rods with tapered free end. It 
is observed that in general the imaginary part of the complex: 
frequency converges while the real part shows little divergence • 
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2 . 

{ 

{ 


The column matrices { 


(e) 


,3 (e) 
n.(ei M~— 

i N 3 

ax 


h 


0 


} >{ 


gje) 

IT. 2 m„ V 2 -4~- 


.3 t 


+ pi 


3T. 'Hi 

i 

4 

Lx 


(x +z) ax 

P 3i 


Keyr.h 


0 


h 
0 

} and 


) 


3 x 3 t 


X=L 


} found in Eqn. (2.14), help 


in introducing natural boundary conditions explicitly. These 
must be included in mass, damping and stiffness matrices. These 
column matrices can not be neglected as suggested by Szabo and 
lee [ 24] • 


3. lor rods with pinned-pinned ends, coupled-mode flutter 
occurs for one particular Telocity only, i.e., at critical velo- 
city far coupled-mode flutter, beyond which ordinary flutter 
occurs. This contradicts the theory presented by Paidoussis and 
Issid [27] that apart from buckling, the conservative systems 
can only be subjected to coupled-mode flutter. It is left for 
the forthcoming researchers to explain. 

4. The results far clamped- clamped rod for old formulation 
are corrected. Early results obtained by classical methods , 
Paidoussis [ 14 3 , show the absence of buckling in the first 
mode. Whereas the results obtained by finite element method 
indicated that buckling will occur in the first mode and the 
critical velocity far buckling has been found to be vezy close 
to the first buckling velocity for conservative systems obtained 
by usirg Euler's method, Paidoussis and Issid [ 27 ] • 
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5. Results for clamped-c lamped rod for new formulation are 
plotted. The behavior of the clampad-clamped rod changes a lot 
by the inclusion of transverse component of tangential frictional 
force. Flutter instability does not occur in this case, only 
buckling instabilitjr occurs in first and third mode . 

6. The finite element method vdth G-alerkin’s approach proves 
to be very versatile and flexible, because it gives a strong 
means to account explicitly for the natural boundary conditions 
and geometrical boundary conditions. 


rod with tapered 
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AEBBHDIX A 

SQUAT I01T OF SMALL LATEBAh LICTIOHS 

The system to be considered for the present study consists 
of a flexible cylindrical rod which has circular cross-section. The 
rod is immersed in a fluid which is flowing with a uniform velocity 
V, parallel to the neutral axis of position of rest of the rod which 
coincides with the x-axis of the co-ordinate system, considered as 
shown in Fig. (2-a). Gross-sectional area of the rod is A, linear 
density (mass per unit length) is m^ ? with flexural rigidity El, 
and total length of the rod is 1. The density of the flowing fluid 
is P f . 

For the formulation, following idealizations are made: 

(i) The rod is assumed to represent a Bernoulli-Euler beam. 

(ii) The fluid flowing past the rod is considered to be incom- 
pressible . 

(iii) Small lateral motions of the rod about its position of rest 

are considered, during which the angle of incidence & and 

remain sufficiently small so that : 

3 x 

(a) no separation occurs in cross-flow, 

(b) the fluid forces on each element of the rod may be 
assumed to be the same as those acting on a correspond- 
ing element of a long straight rod of the same 


63 


cross-sectional area and inclination. 

(iv) The fluid is supposed to be contained by boundaries suffi- 
ciently distant from the rod to have negligible influence 

on its motion. This assumption is made for the system which 
consists of a rod immersed in an infinite fluid and is not 
valid for arod-shell system. 

(v) The virtual mass is equal to P^A for unconfined flow, where 

p is the density of the material of the circular rod; and 
r J 

otherwise C PA, where C is coefficient of virtual mass, 
m r m 

Paidoussis [17] and Chen [ 3 1 . 

(vi) The fluid flow is steady and uniform. 

(vii) The rod is supported at one or both of its extremities, and 
it is assumed that 'die supports are such that they generate , 
at most, only localised disturbances to the idealized flow 
conditions assumed here. 

A . 1 FORCES T0_ EE_ OOKSIIEBED : 

We now consider a small element 5x of the rod, as shown. 

in Fig. (2. a), where T is the tension, Q is the shear force, M is 

2 2 

the bending moment, m^( 3 y/3 t ) &k is the inertia force of the 

rod. F and F_ are the frictional forces per unit length in the 
IT 1 

transverse and longitudinal directions, respectively, and 

the rod gravity. Apart from these forces the rod is subjected to 

fluid pressures and hydrodynamic force which will be discussed in 



detail along with the frictional forces in the next section. 


A. 1.1 FLUID ER EgS UEB : 

Let P and P be the resultants of the steady-state 

px py 

pressure p acting on the outer surface of the element. These 

resultants can be found as follows : Consider the element 5x of 

Pig. (2-a), momentarily frozen and immersed in fluid on all sides 

as in Pig. (2-b); then there will be forces pA and pA + i 

on the two forces of the element in addition with F_ and P • 

px py 

The resultant of these additional forces and of P 3 x and 

px 

P {x is equal to the buoyancy force. Assuming p = p(x) is a 

py 

linear function of x (which is reasonable to assume, as it applies 
to both a hydrostatic pressure distribution and to one modified hy 
skin friction pressure drop), it can be written (see Pig. (2-b)) 


SP = 0 
x 


P 6 x 
" px 


(pA + 6 x) + pA 


- buoyancy force 
= - |£ AS. 


and 


2 

£ P =0 P 6 x - (pA + + ^2 <5 x ) 


6x-(pA +a ^) + 

y ^ py 3x 3x“ 

+ pA f* = o , 

hexice neglecting the higher order terms, one obtains 
- 3? 


px 


.is. a + 2L&1 , 

3 x 3 x 


(A-1) 
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Thus, from Eqn, (A-l) it is obvious that P is equivalent 

P7 

to a tensile force, acting at the cross-section of the rod, of magni- 

3 A 

tude pA. Also, if =0, E =0, hence for a rod of uniform 

3 j - px 

cross-section, the only pressure force acting on the element is equi- 
valent to a tensile force equal to P . 

py 

A. 1.2 HYDRO DYUAmC FORCE : 

The inviscid hydrodynamic force E\ , which is equal to the 
lateral force per unit length acting on the rod, and is given by the 
change in momentum of the lateral flow about the rod, Idghthill i loj . 

The components of fluid velocity in x- and y-directions are, 
respectively, 

V = V cos 9 , 

V = V sin 9 » (A-2) 

y 

where 6 is the angle betv/een the axis of rod and the flow direction. 

Assuming sxnall transverse motions of the rod in either direc- 
tion of its axis, 9 may be approximated by 

© = ( 1/V) I (3 y/ 3 t) + V ( 3y/ 3x) 1 • (a-3) 

Therefore, from Eqns, (A-2) and (/ — 3) one obtains 

y = Y as cos 9=1, 

x 7 

and V = [ ( 3 y/ 3t) + V ( 3 y/ 3 x) ] ■ (A-4) 

y 

Momentum of the lateral flow of the fluid is 

C m m f [ ( 3 y/ 3-t ) + V (3 y/ 3x) ] , 
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hencB ohe rate of change of momentum per unit length becomes 

(3 / 3 t + V 3 / 3 x) [ C m m f ( 3 y/ 3 t + V 3y/3 x)] , (A-5) 

where m f is the virtual mass per unit length of the rod, which is 
given by 



where D is the rod diameter, p „ is the fluid density. C is a . 

X 222 

'Virtual mass” coefficient which is discussed by Chen [ 3] ill detail* 

According to ligh thill C 10 3 , the rate of change of momen- 

tum per unit length gives rise to an equal and opposite lateral force 
on the rod , i . e . , 

F A = C m “f 1 9 / 3t + V 3 /3 z l 2 y. (A-7) 

A. 1.3 FRICTIONAL F ORCES : 

&■ 

Friction forces acting on a rough cylinder set obliquely to 
a fluid stream are taken as proposed by Taylor J 25^ • For rough 
cylinders , he suggested the following expressions for normal and 
tangential drag forces : 


1 

p _ — 

N ~ 2 

2 

P f D V 

2 — 

(C f sin 9 + C^p sin 9), 

(A-8) 

f l = i 

P D Y 2 

C f cos 9, 

(A-9) 


where 0 is the angle between the axis of cylinder and flow direc- 
tion and C„ and C__. are drag coefficients associated with skin 
f JJJr 

friction and pressure, respectively. Hoerner's [8] expressions 



67 


of viscous farces also reduce to Eqns. (A— 8} and (A-9), Paidoussis 

[17] . 

For small oscillations , as mentioned previously) the angle 
of incidence 0 can be approximated by Eqn. (A-3), 

6 = (l/V) t ( 9y/9 t) + V ( 3 y/3x) ] , 
hence Eqns. (A-8) and (A-9) reduce to 

F - — o ETC (AvX , v 
N 2 V^fSt 1 3x' s 

and F l = \ p f DT 2 C f . (A-10) 

A. 1.4 VISCOUS -DA1V1PI AG- FORCE ; 


The force representing the viscous damping effect can be 
expressed as, Paidoussis pl7^j , 



P f D C 


D 9 t 


Therefore, Eqns. (A-10) may be written as 


F. 


If 


P BTC 
f 


f 


(%JL 

K 3 t 




p D c ^ 
f D 9 t 


F 


L 


1 


P •> 3> V 2 C., • 


(A-11) 


A. 2 E QUATI O N OF J IOTIOF : 

The forces and moments acting on a small element, fix , of 
the rod, which has undergone a small lateral displacement y(x,t) have 
been determined and are shorn in Fig. (2. a). The force balances in 
the x~ and y— directions , and a moment balance, Paidoussis [17 3 are 
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H + m r s + *L + (P H + V K - V 


= 0 


(A-12) 


p„ - p. + p + f t ^ + 2. ( T aJL) _ ffl Lz = o, 

I Vi iV Tver T, f\ -V* V v r\ __ / *v> V 7 


3x H 


py L g x 3 x 


3x 


U O 

r 3t 2 


Q 


3 Ji 3 (*!&)_„ (I si_) ) 

3 x 3x 3x 3t 


3x 9 x 


(A-13) 

(A-14) 


where the fact that y and its derivatives are small has been utili- 
zed, and it has been assumed that Euler - Bernoulli beam theory is 
adequate to describe the motions of the rod; moreover, a Kelvin-Voigt 
type of damping in the material of the rod has been assumed to apply* 

Substituting Eqns. (A— 1 ) into Eqns. (A-12) and (A-13), one 

obtains 


yy (T + pA) + (m. g - A) + F y + C? w + Ojf 


9x 


H 


0. 


A 9 X 

... (A-15) 


If + + P&) iff J“ ^2 = 0* 

3t 

(A— 16) 


It - F S - F A + Ft + S 1 Vi - *“■' “r . -2 


How, if yor 0(e), then P^. and are of the same order , 
Paidoussis [ 17 3 , so that 

(f a + P N ) ( 3y/ 3x) 'v 0 (e 2 ), 

neglecting these terms in F.qn. ( A— 15 ) and integrating it from x to 
1, upon combining it with Eqns. (A— 1 5 ) » (A— 1 6 ) and (A-14), one 


obtains 
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yi 




9x 4 9t 


+ El 


^ + p 

3x 4 A 


(T + pA) + m g (L - x) 
jj r 


- |"^ A (L - x) + 


I 

f 


x 


P rtv 



fi + ? u + K « - fi A) 

a z 


3j 

/r 4 *- + m 
ox r 


3 2 

Q — U < 

9t 


(A-17) 


The terms P, E„ and Ih have already been determined in 
A , H 1 J 

sections A. 1.2 and A. 1.3. How the terms (T + pA)^ and will be 
determined considering the case of an isolated rod in parallel flow. 


A. 2-1 CASE Off AH IS0L1TED ROD IN PARAHEL YWV - 

When the cylinder is isolated in a fairly large flow channel 
the pressure drop due to flow may be considered to be negligible, 
and 


3P. 

9x 


P f g 


(A- 18) 


Considering the term (T + pA)^ s if the rod is free to slide 
axially (or is free completely), at x = L, one obtains T^ = - P^ 
(since pA arises from the pressure on the sides of the rod), where 
is the base pressure? hence (T + pA)^ = (p^ - P fe ) A. Since 
P b < P L generally, this represents, more or less, "form” drag, more 
accurately referred to as base drag, which may be expressed as 

(l + pl) T =l p. dVc.. (A -19) 
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If the total length of the rod is fixed between fixed supports s 

1 2 

at x = 1 there will be a compressive load - p, D V (l/2) 

-(m g - A) where the existence of the term m g(l/2) depends 
X* Q X C. X* 

on w/hether the length was fixed while the cylinde r was horizontal 
or not. This fact can be explained as follows: 

Eqn. (A-15) can be rewritten after rearranging and neglect- 
ing the higher order term as 

Tx + P A ) = “ ( m r S ~ U ~ \ ' (A-20) 

Integrating Eqn. (A-20) between the limits x to L one obtains 
tensile force distribution along the whole length of the rod as 

(I + pA) x = (T + pA)^ + ( m r S ** A ) - x) + - x) - 

(A-21) 

Now Hie stress - strain relation gives the strain in the rod 
due to tensile forces 

£ ■ a ♦ *>*• (A ~ 22) 

where is change in the length of the rod. 

Substituting Eqn. (a- 22) in Eqn. (A-2l) and integrating the 
resultant expression between 0 to L, one obtains 


?0 n= k f O [ (T + ^>1 + °"r g "If A) (I ' x) 


+ L ( I - x) ] dx 

JJ 
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°r 0 = ~ t (T + pA) L I + (m, g - A) (L 2 - + F^L 2 - f- ) 1 • 

U-23) 

since fixed length requirement implies that must be equal to aero. 
Solving Eqn. (A-23), one gets 
(T + pA) L = - P l (1/2) - (m r g - g A) (l/2) 

= -[ PT P f D ?2 C f + m r S ( L /2) . (A-24) 

Moreover 5 there will also be a compressive load due to radial 
contraction, for a thin tubular rod this is equal to - 2v ( pA ). 

And there may also be an externally imposed uniform axial tension 

T. 

Therefore for this case 

(T + pA) L = -[ \ P f D Y 2 C f + m r g - g A ] (l/2) 

+ (l-2 v)(pA) + T . (A-25) 

CombiniJig Eqns. (A-19) and (A-25), and substituting into Eqn. (A-17) 
along with Eqns. (A-18), (A-7) and (A- 1 1 ) , one obtains the equation 


of small lateral motions 
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V 


? fl- 
ax 4 at 


+ EI + C m 
3x 


, a _ a %2 
m i C 5t + T 3^ y 


_{«[¥ + (1 - 2 V ) ( pi )) +[■! P f W 2 C f + ( m r - 

[ (1 -IS )L-x) *\ P £ B 2 '.' 2 (1 -S) \}^i 

( 1 X , v -- p D C CL |^r + (m - m ) g 
l 3^+ v 3 X J + 2 f D at r X ax 


m f ) g ] 

+ | Pf D V C f 

+ m - 3 “! = Oj 
r 3t 2 

(A- 26 ) 


where 6 = 0 signifies that the downstream end is free to slide 
socially, and 6=1 if the supports do not allow net axial tension. 
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APPENDIX B 

E QUATION OP MO TION (OED E O MUIATIOl) AM) 15 'S FINITE EIEiENT 

ANALYSIS 


Neglecting gravity, pressurization effects, inrernal damping 

and the y-component of tangential frictional force F T , i.e., F (~£) 

L Lax 

Eqn. (A~26) becomes 


“ -H- + “V (5 * + V A) 2 y - { S I 




m f at 




2 p f 


1 p 

2 f f 


BfC f l [ (1 -6 /2) 1 - x3 } 

3 x 

MC f ft + v l$ + i p f D?2 c f li 


+ m = 0 

r 3t 2 


(b-1) 


Eqn. (B-l) represents the equation used by Paidoussisl 14 3 
which has been referred as "Old formulation" in the present work. 

It docs not take into account the frictional forces correctly. Com- 
paring this equation with Eqn. (A-26) , one observes that in new for- 

a 2 y 

mulation the term obtained by combining the coefficients of x ? 
and is F_ (x + P3 (neglecting gravity) and in the old 

O X Jj n ^ O X O 

ax 3 T 3 Y 

formulation the aboYe term becomes “2 * ^ " fciie neW 

3 x 


formulation, integration by parts of the above ternn is straight for- 
ward, sec Eqn. (2.15). Integration by parts of the corresponding 


term for old formulation was not so straight forward. This diffi- 
culty was overcome by splitting this term and then integrating by 
parts as follows: 
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P ( x + Pr) = i 1 , ( x ■Mj -4 P) + B pL . 
L 8x 2 9x L 3x 2 9x L3x 


(B-2) 


Integration of all the terms of old formulation is given in 

9 y 

Chapter 2, except of the term B^^y . Applying finite element based 
Galerkin’s method, this term in Eqn, (B-2) turns out to be after 
integrating by parts 


h a J e ' 

f N. 

J i Box 


dx 


B i F i y 


(e) 


h dl. 


0 


dx 


B^ y^ e ^dx 


(B-3) 


Substituting Eqn. (2.25) into Eqn. (B-3), one obtains for the element, 


nn (e) r , ne 

- [k4] {y> , 


^ 


where [ k4] 

h 

= h > 

{n’} L 

I J dx 



0 




“-1/2 

h/10 

1/2 

-h/10 

s 3? r 

-h/1 0 

0 

h/10 

-h 2 /60 

L 

-1/2 

-h/10 

1/2 

h/10 


h/10 

h 2 /60 

-h/10 

0 


(B-4) 


(B-5) 
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and the "boundary condition oo be considered becomes for the whole 
domain 


- 1 /2 P f D V 2 C f y (e) 
0 
0 




s, 

/ 


0 


+ 1/2 P f D v 2 C f 


m+1 


(3-6) 


It should be noted that the values - 1/2 P, D V 2 and 1/2 D V 2 n 

XX if 

are summed up with ^ and ^2m+1 2m+1 e ^- emen ^ s °- "the assembled 
stiffness matrix, respectiveljr. 
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APPENDIX C 

B0U3HXARY COiflDITIOHS FOB A CANTIIEV^^ D HOD 'ffiTmH TASSBED EREE END 


If the downstream end is free-, terminating in a tapering endj 
the cross-sectional area of 'which varies smoothly from A to zero in 
a dis l.-ance 1 ( << l) , y and the lateral velocity v may he consi- 
dered constant for the length 1, Paidoussis [14] . This requirement 
allows the forces acting on the tapered end to he lumped and consi- 
dered in appropriate boundary conditions, Hawthorne [ 7] . In the 
transverse direction equilibrium equation over the tapered end, say 
for 1-1 x < 1 , is 

f cte - f f ( sf + 7 -~) E m (x) v ] dx 


1-1 


1-1 


1-1 


m (x) dx = 0 

T rv d 

3t 


(o-D 


The parameter f , which is equal to unity for slender 

0 _ 

budy, inviseid flow theory, has been introduced because the theore- 
tical lateral force at the free cud may not be fully realised as a 
result of (a) the latei’al flow not being fully two-dimensional, since 
the fluid has opportunity to pass around rather than over the tapered 
end [ Munk , M.M. , 1924 MACA Rep. Ho. 184 ], and (b) boundary-layer 
effects, Hawthorne [7 3 * Accordingly f Q will normally be less 
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than unity, Paidoussis [ 14] . 


Writing m„(x) 


A(z) and m^(x) = A(z) and, since 


y and v are constant over the tapered end, Eqn. (o-l) leads to 


1-1 


3x 


dx 


' 1 p i T T ! si - ( p r + 1 p ± ) ri / 


L-l' 

A (x) dx = 0 


which yields 
El 


where 


3 3 y 


>4 


3 


X 


3 + PI tKi + f = m f T ( a"t + T 5 - 


(m + f m ) xe = 0 , 

v c f ^ 


a t 


xe 


1 ^ 

r J A(x) dx . 


L-l 


'l-l 

( 0 - 2 ) 


(0-3) 


(0-4) 


Eurthormore , it may be assumed that there is no bending 
maaent at the freo end, thus 


3 2 y 

El -H 
3x'‘ 


pi - 


3"V 


= 0 » 


(0-5) 


3 x ' 3 t 
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APPENDIX D 


THB COMPUTER PROGRAM 

The computer programme was prepared for assembling the 
system matrices, substituting the boundary conditions, in-vert ing 
the assembled stiffness matrix, forming the dynamical matrix and 
finally for finding out the complex frequencies of the system? 

lhe computer programme was written in Fortran IV and run 
on an IBM 7044 computer with approximately 32K words of core storage? 
The programme is quite general and built-in subroutines MATINV for 
matrix inversion and ZINDLA for finding out complex frequencies were 
called directly. 

Fortran listing of the program is given in Appendix 
D.1 DESC RIPTION OF PRQSRAM ; 

The MAUI programme reads and pi’ints various input data and 
output results . Other subroutines calculate the matrices, assemble 
them and after applying boundary conditions form the dynamical matrix? 
The complex eigenvalues are calculated by calling built-in subroutine 
a INDIA. 

S ubroutine Na me __ Description 


STIFF 


Called by MAIN, calculates the 
various elemental stiffness matrices. 
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Subroutine 


MASS 


DAMP 


ASMBL 


BNDRY 


MATMPY 


Name 


Description 

total stiffness matrix, assembles 
bnem, introduces boundary conditions 
by calling BIIDRY and calculates 
inverse by calling MATIffV. 

Called Dy MAH, calculates the ele- 
mental mass matrix, assembles and 
introduces boundary conditions by 
calling ASMBL and BNDRY subroutines. 
Called by MAIN , calculates the ele- 
mental damping matrices, sums them 
up, assembles them and introduces 
boundary conditions by calling ASMBL 
and DNDRY subroutines in turn. 

Called by MASS and DAMP, assembles 
the mass and damping matrices in 
turn. 

Called by STIPP, MASS and DAMP, 
introduces boundary conditions to the 
assembled stiffness, mass and damping 
matrices and separates the required 
matrix:. 

Called by MAIN, forms the dynamical 


matrix, 
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Subr outine Same Description 

MATINV Called by SHE'D, computes matrix 

inverse . 

BIGENP Called by MAIN, determines eigenvalues 

and eigenvectors of a general real 
matrix. The binary version of this 
programme is stored on 1301 disk file, 
which is named Z I1ID LA (Dor further 
details refer to communication of 
ACM journal, Vol 2, 1968, pp. 820). 


LIST QD PRINCIPLE VARIABLES 

FORTRAN Dap lenient at ion Def initi on. 

P rogram Sym bol 


MIN PROGRAM 

AL 

AMI 

AMP 

AMP 

AREA 

AV 

BETA 

CBAR, CRASH 
CM 


length of the rod. 

Moment of inertia of the rod. 
Mass of fluid per unit length. 
Mass of rod per unit length 
Cross-Sectional area of the rod 
Actual velocity of the fluid 
Dimensionless mass of fluid 
Driction drag coefficients. 
Coefficient of virtual mass. 



FORTRAN Impleme nta tion. 
Program Symbol 


Definition 
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ilA.Il PROGRAM 
DO, DI 

E 

EVBl(l) 

EVRl(l) 

FAC 

FI 

I 

1MI1 

NJUAZ 

ISP 

OlfiGAI, OM3GA2, 01BGA3 

ROW 

RCWP 

VNON 

SVIOH 

SUBROUTINE STIFF 
AXP 


Outer and inner diameters of the rod, 
respectively. 

Modulus of elasticity. 

Imaginary part of dimensionless 
e igenvalue . 

Real part of dimensionless eigenvalue. 
Dividing factor which gives dimen- 
sionless e igenvalu.es , 

Tangential drag force, 
lumber of elements. 

Minimum number of elements. 

Maximum number of elements. 

Total number of velocity points. 

Lowest three exact natural frequencies. 
Density of fluid. 

Density of rod. 

Dimensionless velocity. 

Dimensionless velocity data. 

Axial distance from upstream end to 
node p, see Fig. (l- c). 
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Program Symbol 
SUBROUTINE STEPP 
ELME 
PACTOR 

mu 

N 

HI 

P 

GSB 

GSM 

S 

51 

52 

53 

54 
S6 
X 

z 

SUBROUTINE MASS 

AM 

GMB 


Definition 

length of one element, 
largest element of GSM 
Factor multiplied to elemental 
stiffness matrix XNDHD (4 S 4) 

Size of the assembled matrix = 2 x (i+l) 
Size of the separated matrix 
Factor multiplied to elemental beam 
stiffness matrix 

Assembled stiffness matrix (U,n) 
Separated stiffness matrix (Hl,Hl)* 
Total elemental stiffness matrix (4 ? 4)» 
Elemental beam stiffness matrix (4,4). 
Elemental stiffness matrix MFVSQ (4j4). 
Elemental stiffness matrix DY/DX (4,4) • 
Elemental stiffness matrix XMDHD (4 S 4)* 
SI + S3 - S4 

Factor multiplied to matrix S3* 

Factor multiplied to matrix S2 

Elemental mass matrix (4j4)» 

Assembled mass matrix (H,H). 
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Program Symbol 


SUBROUTINE MASS 

gmm 

H 

P 

SUBROUTINE DAMP 

alpha 

C 

Cl 

C2 

C5 


GCB 


COM 

N 


P 

PIRKT. 

SUJiROU'iM .'L ASi'.LL 

AA 

ASM 

N 


Definition 


Separated mass matrix (lTl,Nl). 

Size of assembled mass matrix=2x(l+l) 
Pactor multiplied to matrix AM, 

Pactor multiplied to matrix C2 
Elemental damping matrix AMFAV (4,4) • 
Elemental damping matrix DY/DT (4,4), 
Elemental internal damping matrix (4,4), 
Total elemental damping matrix (4,4), 
Assembled damping matrix (N,1T). 

Separated damping matrix (Fl,Fl). 

Size of assembled damping matrix 
= 2 x (l + 1). 

Pactor multiplied to matrix C. 

Pactor multiplied to matrix 01* 


Elemental matrix (4,4). 

Assembled matrix (N,N). 

Total number of degree of freedom 
of the system (size of ASM)=2 x(l+l). 
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SUBROUTINE BUDRY 

ASM 

BSM 

N 

N1 

SUBROUTINE lUATMFY 

D 

3)1 

DZ 

mm 

R1 

NB 

SUBROUTINE MATIHV 
A 

B 

M 


UBTKKK! 


SUBUOmu.'iE Fi.UlEili 1 
A 


Definition 

Assembled matrix (H,H) . 

Separated matrix (ill, Hi) 

Dimension of ASM = 2 x (i-fl). 

Dimension of BSM. 

Unitary matrix I , (Hl,Hl) 

I'k " 1 m] s (II, HI) 

[k” 1 cj , (HI, HI) 

Dynamical matrix (HB, HB). 

Dimension of Qcj , jjs] , [c"j , 

M md [°J 

Dimension of dynamical matrix = 2xHl 

Coefficient matrix of size Hi to be 
invorsed. 

Constant vector of size HI 

Humber of constant vectors 

(if !!s 0) only inverse is computed). 

Value of the determinant 


Matrix of size Hi x Hi for which 


eigenvalues and eigenvectors are 


determined* 
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Progra m Symb o 1 
SUBROUTINE EIGENP 
EVI (I ) 

EYR (i) 

MS 

T 

VECI (J,i) 

YECR (j,l) 


Imaginary part of eigenvalue 
Real part of eigenvalue 
Dimension of dynamical matrix given 
in the main programme ( 1 'JM^.N) 
Dumber of binary digits in the 
mantissa of a single precision 
floating point number = 27,0 
Imaginary part of eigenvector 
Real part of eigenvector 
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$ 3 ? '$ $ $ $ $ 4% % 


. DATA CM,.C£!.T«/l.-j 
*'■• DATA. On/’ . / 

S £ e GAFA~1,av : 
s : T = (GAVAvt=*Ay I ) / ( { A L 
%. P R I' : 7 154, C- A i v A , G 

!|®54 . FORMAT (7 , "5X,>GAM A=* 

► CGAMA*0. :» 

lie . TBAR=<CGAMA*'.*Af11)/(AL**2> 

Pi. . PR INT 155 , CCA 1“ A , T 6AP 

§St55 FORMAT! 5X,*CCAMA=*, :; 16,6, 1CX, *T3AP.=* , €16 ,8 } 

S $i $ $ if if if if jjc ^ ^ $ if if ix if j^c >e ^ 3^5 i&if if if, 3^5 3^s if~f ^j* ?$c sjc 3$t ?$c <{jc )$c 3$c 5^c 3^£ if, 3§c ?*£ if 5^ s$s if, 3$c i£c 

Is ' SATA FT A j Mi / «n, '•.::/ 

lie: PAEAi=f cTA*E*AMI )/UL**2) 

*'” INT156, ETA«NU, PA EAR 

R M A T 15 X , * E TA=* ,F 1 6. 8, 5X, *il=*, 016. 6,5X ,*PA3AR=*,E!6.8) 

*#****#•*•#**** ******* ********:?******************************* 

MT*A1/FLCAT( I) 

********** ****** * *********** **************** ********** ********* 
[6*MI >/!EL'MT**2) 

ti2J*6^eLMT*P 
1,3)=-Sil, 15 
l,4')*S(lt 2) 

2»2)=4.*{ELMT**2)*P 
■2«3)*-S! 1,2) 

2,4) = 2.*!ELMT**2)*P- 

tltKIhH, 

AfA)*-$('2»2> 


AFF) 5 


1, JL ) 

* J JJ *GE* 3) GO TO 78 

lliAH ST I F FNE S S m 1RI X IS*) 
;ill®l)t JL = 1,4),II = 1,4) 


*************************************************** 







7 (■ / » 5- X . : 
: ? » ( ( s * t : 

! / « *i X 5 4 


: FN ?S MATRIX 

j K J » J K = 1 j i ) t ik = : 

i X ) 


»*« *5* t,.f 4 v^* Vrf «* ' 

*v ^ v *4'' nr -v* «v 


fi^S s*? & $$$$$$$:$$; 


■w 



S I (2 
r - 1 ( 2 
S 1 { 3 
SI (3 


SI (4 1 1 ) = S 111*2) 
i S i '(4 , 2 )=~S 1 C ?,4) 

: ■ S l (4» 3 )=-S 1 f 1,2) 

.$114,41*0.? ' 

Iff f .KE.2.CF.JJJ.GE.3) GO TO £1 

K -PRINT ISA 

814 '.FORMAT { /t 3X* ^STIFFNESS MATRIX DY/DX IS* ) 

® ' RR1NT23? ( ( SI < IK » JK), JK = 1,4),IK*1,4) 

B FORMAT t/,4X,4 c .2-‘,F) 

pisiiifiiwe" 

Sp $ sf: S^c sp 3^$ 4$ $$$$$$$ A $$$$$$**$$££ $4$ 

oiTi { It=|.,4 ' . . . . ' 

., K p|||^aU=S6 { IL, JL )+{ S ( IL, JU + S1 (lit JL )-S4 l IL, JL) J 

|lfc||l*4^Fi‘UV**2)->{DELTA*ITBA'R+.(l. )-2.0*NU }*{ PABAR ) )+W* { { I 
c!i||t|)! + ^^^ CU ’.-.I* 5*DELTA)*AL)+f0.5*R0WF*(D0**2)*(AV** 


||LMT**2}*Z 



21 

7° 


41 ' 

; r-4 

.24 

82 ■' 


C 7 f 

: , 4 5 =■ 

2 f 

0 , ; } a 

. ; “ i 

3 ; : : } = 

~ ; } 

^ f l ^ 


JK 

: ;• ( 

K , ! K 

i f { 

t •» - 

i e ' - «• 

■ . o '; A 

/i ' : 


o, > / 

■ 1 ’ ' ■ ! 

'r 0-; IV 

: ; ( / 
7 T » f ? - 


)*: 


( ,2) 

: » 2 ) 


■ j. ? ■ 

: i j 


f IK, JK 


f ■ 

f ( 


f y - b 7 4- r T~ j‘: •■. S b A T 7 

*( IK, JK) - JK=l ? 4), IK 

c- 5 } . 


TO 79 
X WFVS< 


7 , * 


r 4 ) 


, S $ ) 


i • ’ ~ i. ? 

J K = 1 * L 


c -s f 

* 

% ■ * 

7 C / 


, J 

'ill 


St { IK, JK }* S2 C IK, JK } 
CF.JJj.O *. *1 r.t 


3) GO TO S2 
F FT ESS MATRIX 


> 


e: -:”? v A' ,, (/,?K f *TCTAl 

3 Ri j 'T24 , { (S3 ( IK , JK ) , JK = I , 4 } , I K =1 , 4 ) 

" JlM.57( /* 4X, 4~2 .H 

0 0 V 7 I 3 U '• 
r. f - 3 | i 

J J = 1 1 

1 13.* T 1+3 
J J 2 = 2 J+3 

TO 1 IK=IT ,11 3 

00 122 JK=JJ,JJ 2 
L = IK- 1 1 + 1. 

LJ=JK-JJ+1 

GS E { IK, JK)*GSBIIK, JK)+S3CL,LJ> 

OOMIMJ? 

CALL BNCRY { G'SB, GSM ) 

.... .. .. -•-*--•••*-> . -< ¥**»***■*.• 

: a= iz*ii+ n- 2 ) 

F.A-CTCP=ABSCGSM( 1,1) ) 

SO 9 IK=3,N 
DO 9 JK=1,L 

1 IFIAfiSI GSM (IK* JK ) ) .01* f- AC TOR ) FACTOR =A 8 S(GSM(IK»JK) ) 

’ CGNTINOS- 

; I F f I * K .• 2 « G R0, J J J « G E3 3 ) GG TO 83 
J>R|NT37, FACT Oft 
: P G m A T ( 5 X * * P £ C T C F = * f c 1 1 . 8 ) 

CONTINUE 
300./ 13 : I K-1, P 
06/13. J K = 1 , Li- 

’OSffTIK, JK 1=6 SM ( I K, JK ) /FACTOR 

wxsoe ■' 

m i#£*2«CR.JJJ*GE.3) GO 10-84 




3X^ * ASSEMBLED STIFFNESS/FACTOR IS*) 


mmm f- \ / * ■* a i i r r i\ c j o / r wu i u 
(GSM IK * J K ) , J K = 1 ,N) , IK = 1,N) 

; iMT-f/:,:5X,fcE15.8) 

IfTfWS' '/ 

lit WS- ' SUBROUTINE M A T I N V . 

v ■ 


ifill* 

//•I 


tiafiSsUL'Si'SMB’M 


ERSE OF GSM IS*) 

, JK 5 * J K=1 * N 3 * IK — 1, N) 
• 8) 



I* K J ) / FACTOR,; 



■* J J J 




t '-2 l . ~ } 





s 5 






-C0>*Cf-/ei3/A'Fj : i.MF>FL 
. cc^r/eH/R 
'. C C>MC‘. / : ?1 5 / ft CWP> PCWf j AF. t A 
I • , ■CC*'MCN/::i! 7 /I , JJJ 
lift DIMENSION GVH 2 X t 2 S) 

I CIS' AAtK f 10 >,A$M< 2 C, 2 ' ),£SM {2 

/ F L C A T ( I ) 

s$.#*3jt*a5!.jfrst ^ afcjJ' ***£** sJsjSs^ki^S: :£ 

4 AMn*ELM 7 >/ 420.0 
15 fe»*P 
22 **E LMT* P 
54. *P 

~ 13 .*nLMT*P 

4 |||RtMT** 2 J*P 

~ 3 .*( ELMT ** 2 ) *P 

!gE. 3 ) Gfj T n €7 

/? 5 JU*«ASS MATRIX IS*} 

|TI|| I I » JU,JL = 1 » 4 ),II = 1 , 4 ) 


■ # 4 *** ************************************* ******** 
&#■****** ************ **.*** At* ******* ********* ******* 

it^OAM-p TNG .Matrix amfav. 

TAL DAMPING MATRIX DY/DT. 

TAL INTERNAL CAMPING MATRIX. 

ELEMENTAL DAMPING MATRIX. 

BLED DAMPING MATRIX. 

AT EC DAMPING MATRIX. 
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I ; r t 


( I ) 




«A» ’-V «J 

•'* *•<* - - *1- *.!« AT 


¥ 5 I : V : 


1 "V* V V *(* Sk 




c- 
( ': 

I ". 

ll 


r, 

r 

I-r 


• r, 3 — ... 
. 4 ) =- 

t 1 ) “~ 
f it ) - 

,0=3 
»<•■’ = •* 

, * ) =~ 
, ; ) =.-, 

» ~ ) = -* - 

,4 ) = 8 

s I 5 = t 

, n = ( 


- : l < 

i >» 


OF 


( 4 

( A * 3 ) = -(. 
(4,4 } 

” ( I « : 
i 


* r u-7*e 
ftn m*p 
LFT ft$2 )ft? 
.ft&LNT*£. 

p* jjj*o 


3 5 GC 


88 


lumping matrix is*) 

JK = I,4), IK = 1»4) 


!•' AVftCBAF 5+COASH) /42 u 
K T 


e £ © f. 

f K 4 ’■ i A X 

FO.P>AT{ /,F*,*FI.OT 
« :, F ,((£.{ IK, JK3,. 

F Q.P ?■:' A T ( / , 4 X , 4 F 2 G o 8 ) 

C:Q K7 I rU-3 

4 * * * ftftftft ft ft ft ft ft ft ft ft ft ft ft if ft ft ft * ft ft * ft ^ ft ft ft ft ft* : 

P IN- K I a ( { 3 Cfc F* ( DC / .? * ) * AV : ' ~ ~ 

c. X ( l , 1 ) a 5 h t , * 3 L H T* P I !' . K I 

Cl ! I ,25 = 22 »* ( CLFTftft ) ftp INK I 

G 1 { 1 J ^ ) = ''A '* CM HtT-J.PT’ V T 

C1C 1,4 )=-o‘ 

Cl (2,2) 

C 1(2,3 5^ ........ 

0 ( 2 ,4 ) =-3,* ( FLFTft 
0(3,35=0(1,1) 

C l ( 3 ,4 5 -~C K 1 » 2 3 
Cl (4 ,4 5 =0(2, 2 5 

ft.ftft ft ftftftftftft ftft ftftft ftftft ftSftftftftftftftftftftftftftftftftftftftftftftft^ftftft.ftftft.ft.ftjJ,^ 

DO 4 I K= 1 , 4 
0024 JK=1, 4 


*-ELPT*P10KI 
="• O * * ( F L FT ft * 2 } ft p I C K I 
=4*ft { F L Mftft ^ 


-C l ( ! ,4) 


■ ft 


) "-PICK I 
5 ft P I N K I 


m 4. JK'= 1 , 4 

i,*.-.--.vltaK,lK)aCl(lK 
... IF(UHE.2,CP " 
■ PRINT 111 
' FORMAT ( /,5X, 
PRINT 31 


,j„,JK) 

• JJJ.GE. 31- 


GO TO £9 


, * : 


•r.t*ws m n / f d a * * 5 r; C -CN D ' D A M P IN 6 MATRIX 15^) 

2 P m « T 3 1 * n C i. U K , J K ) , J K s 1 » 4 ) > I K * 1 , 4 ) 

O' . F OPM AT ( / » 4X , 4E 20 » £ ) 

ftftftft ftftftftftftftftftftftft ftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftcftftftftftftftftftft ftft 

f A N P + A F F 3 / A M I) 5 ft A L P H Aft ( AL ftft 2 3 
H 1350 , AL PHA , FU 

tPHA=ft f £!£• 8, 6X, ftHU=ft, E16.8 3 
ililflELNTft^l 


■: i 

O: (1, 3 3=6*47 LHTftp 

■V : CTT liBl^-CZC 1,1) . 

T- | — .f* 3>.|: | y 2 ) 

" |||44ift.( £LMTft*:21ftP 


non 


\ : * .. L -- 


g \ ... • i J 

C {• ; , 1 ) 

•: f ; „ v 


(Tit Jl ) 


: * * « V * *1 ; ' # ~ *> % * * * - ## * # « # 5$ t * *- % # ± ^ * # * # $ * $ * # * A £ £ 2& # J* .: 


■: { 1 L 1 J L 5 : 

f * * , . I * \ ; 


' t 


' r • * 
r :t. *■"' 


r 2 f IL, JL )+ (C CIL , JL. ) + CKIL, JL )+C2 { It, JL }} 
••jJ J*X ':) GO 7' ?: 


f / ? " T C T 2 1 Y & p p I f\ p p 7 R 1 >' I $ * ) 
■ - * n ~ ~ ■ - 

• ?■ 5 "' 


;; * t ( c 3 ( IL* J L ) , JL = 1 , 4 ) , I L = 1 , 4 ) 

•• : ( / * 3., ;/ * tv 


: 'X'XL G C E ) 

: - 2 Y { G C R f G C ’ ) 


. 2 . 


i 


■XX 2 

ii 


:- u r : 


,:■ "■ r ? • 


T :i * * A * * * - * *$$* 4 4 **** * * * *ft £. A A * * A * * * * A A * * * A A # * * * A ^ # A * * * 

X A XXL (A A, ASM 

'' ~ 5! - J n- tfXXS* -? S? A * t- * * Jjt # * aje sjc ft * * ftft * * A if AAAAAAAAAAAAAAAAAAAAAAAAA # 

■ ’ :v XE!\T4L MATRIX. 

R I X» 


7/ 


L'FHPl 

» JJ.J 


) 


> f 


, : ; 4 ; A 

» A S X { 2 


2-R),R5M{ 20,20) 


;V f 


;p J '* 

;■• p ft. , f # p V j 

: ) T f v - :v :; t 

( T + 

0 ' : T T 

::kj i Jl-ifP 
^5K{ II ? JL) = 

OC 2 :■< - . ? I 

1 l ~ 9 . ❖ K — * 

) L = I I 
! 13=11 + 3 
J 13= JL+3 
90 2 !K=I] 

3G 2 JK---JL, 

L = IK- I T + - 
L J = JK--JL+ 

ASM IK, JK )~ 

I F ( I « M * 2 » 

. PR 1X7 1 1; 

103' FORMAT ( / , - X , : ^ASSS y BL “D MATRIX AFT-R B 
PR INTI J. 3 , ( (A-SM IK, JK), JK = 1,N), IK = 1,N ) 
1S4 FORMAT (/, fX, 671S.P) 

79 COM INUF 

'■ P'- p T pOfo ■ 

" : EKZ 


SUERGUT IN T: I’MRY (ASM, 8SM ) 

mn IS. T H 2 ASSEMBLED MATRIX* 

SM IS . THE SEPE BATED MATRIX. 
I3MMW817/I,JJJ „ ^ ' 

-pgiSieN-' ASM ( 20 , 2 C ) , 8 SM ( 2 C, 20 ) 

ARY CCMIITICNS FOR PINNitC-PIMMED ROD# 
*(! + !) 


i! 


SM I K , J K ) + A A I L , L J ) 

■ rc *JJJ.SE,3) GO TO 79 


C. IS*) 



mifh 

Wf*QW-. 






I L = 

jK=2,n 

j L - j ; + 

lit JJL )=AS ; N ( IK, Jk ) 

, * -•» t J. - v, .■ 

- ■- f J . ■ J 
'•= ? =!' [ 1+ ) 


! F ( I * 

.pp. i::t 


f *J J J < 





* ^ $ 5?C $ # 3{C # $ '* 4 .■# $ J* # * 3$. . $ ** ** * * $ * ^ ^ ^ * $ * 3$ 4 * * * 


0 IS the IM TAPY f'AYRIX. 

01 = 1KVlHS :: - C F ,GS w -*GXV. 
f >2 = IKV'.^?t. 1 CF GSM*GC-M# 

!>MM IS TFF CVNAMICAt MATRIX 
C'CMNCN/S] 7/ T i J J J 
COMMON/*? /GSM 2 r .:,/G} 
C0MM€K784/G^V.{2v»2C ) 
CQMMCN/66/GCMt 2 9,20 3 
COMMON. /B 1 8 /l M M ( 3 C * 3 €■ ) 

'COMMON /B9/0(7n»2S) '' 

COM MCN / 6 1 9 / C I ( 2 w » ? 0 ) 
COMMC?\/£8/C:'{2C,2f ; ) 

^ UNITARY MATRIX D* 

>‘C ( | 

||i|||l:BF) 0( IMvJH)=l.C 



: , JL )+ (“• GSP ill ,K)*GKM K, Jl ) ) 


L , J ) + (-GS S ( I L , L ) * GC II ( L f J L } ) 

CAL FI IT I>; r f ' ’•< c 
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^(r.j;')TT! IK, Jl ) + D I ( II, JL 3 
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J }.- f V 


■; » [ I’< , -JK 

)« 

IT( IK j 

JK )+! 
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3 ‘J J J , •* 
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